classdef ParametricBlochEquation < Algorithm.BlochEquation.BlochEquation
    %PARAMETRICBLOCHEQUATION Summary of this class goes here
    %   Detailed explanation goes here
    
    properties
        x
        kc
        
        pMax
        nMax
        nT

        
    end
    
    methods
        function obj = ParametricBlochEquation(T1, T2, varargin)
            %PARAMETRICBLOCHEQUATION Construct an instance of this class
            %   Detailed explanation goes here
            
            p = inputParser;
            p.KeepUnmatched = true;
            
            p.addParameter('ratio_x', 1.0, @isnumeric);% ratio of 'control amplitude' to 'control frequency'
            p.addParameter('ratio_kc', 1.0, @isnumeric);% ratio of 'control freq' to 'omegaZ'
            p.addParameter('pumping_rate', 10e3, @isnumeric);
            p.addParameter('tMax', 100.0, @isnumeric);
            p.parse(varargin{:});

            Rmax = p.Results.pumping_rate;
            tmax = p.Results.tMax/(Rmax + 1.0/T1 + 1.0/T2);
            
            obj@Algorithm.BlochEquation.BlochEquation(T1, T2, varargin{:}, ...
                                    'steadyOpt', 'manual', ...
                                    'tRange', [0 tmax], ...
                                    'steady',    [0 0 -Rmax].' );
            
            obj.x = p.Results.ratio_x;
            obj.kc = p.Results.ratio_kc;
            obj.R = Rmax;
        end
        
        function sol = solveAnalytically(obj, varargin)
            p = inputParser;
            p.KeepUnmatched = true;
            
            tTest=@(x) x >= 1;
            pTest=@(x) x >= 1;
            nTest=@(x) x >= 0;
            chebTest = @(x) isnumeric(x) || isa(x, 'chebfun');
            
            p.addParameter('nT', 100, tTest);
            p.addParameter('pmax', 20, pTest);
            p.addParameter('nmax', 20, nTest);                        
            p.addParameter('omegaX',  obj.omegaX,      chebTest);
            p.addParameter('omegaY',  obj.omegaY,      chebTest);
            p.addParameter('omegaZ',  obj.omegaZ,      chebTest);
            p.parse(varargin{:});
            
            obj.nT = p.Results.nT;
            obj.pMax = p.Results.pmax;
            obj.nMax = p.Results.nmax;
            
            obj.set_Omega('omegaX', p.Results.omegaX, ... 
                        'omegaY', p.Results.omegaY, ...
                        'omegaZ', p.Results.omegaZ);
            
            omegaC = obj.kc * obj.omegaZ; Tc = 2*pi/omegaC; % control freq. & period
            
            gamma = obj.R + 1/obj.T1 + 1/obj.T2; % rad/s, energy unit
            tau = 1.0 / gamma; % s, time unit

            bx=obj.omegaX/gamma; by=obj.omegaY/gamma; bz=obj.omegaZ/gamma; zc=omegaC/gamma; 
            
            tList=chebpts(obj.nT, [0, Tc]); 
            obj.t = chebfun('t', [0, Tc]);

            [mx, my, ~, ~] = obj.mxmy(bx, by, bz, tList/tau, obj.x, zc, 'pmax', obj.pMax, 'nmax', obj.nMax);
            
            obj.solution = cell(1, 3);
            obj.solution{1} = chebfun(mx.', [0, Tc]);
            obj.solution{2} = chebfun(my.', [0, Tc]);
            obj.solution{3} = chebfun(1, [0, Tc]);
%             obj.solution.mx_p = mx_p;
%             obj.solution.my_p = my_p;
            sol = obj.getSolution();
        end
        
        function sol = solveNumerically(obj, varargin)   
            p = inputParser;
            p.KeepUnmatched = true;
            
            tTest=@(x) x >= 1;
            p.addParameter('nT', 100, tTest);
            p.parse(varargin{:});
            
            obj.nT = p.Results.nT;
            
            obj.t = chebfun('t', obj.domain);
            omegaC = obj.kc * obj.omegaZ;
            obj.set_Omega('omegaZ', -obj.cosine_pump(obj.x * omegaC, omegaC/2/pi, obj.omegaZ));
            obj.set_pumping( obj.R ).solve();
            
            period = 2*pi/omegaC;
            nP = floor(obj.domain(2)/period);
            range = [ (nP-1)*period, nP*period ];
            
            obj.getLastPeriod(range);
            sol = obj.getSolution();
        end
        
        function obj = getLastPeriod(obj, range)
            Tc=range(2) - range(1);
            tList=chebpts(obj.nT, range); 
            obj.t = chebfun('t', [0, Tc]);
            mx = obj.solution{1};
            my = obj.solution{2};
            mz = obj.solution{3};
            
            obj.solution = cell(1, 3);
            obj.solution{1} = chebfun(mx(tList), [0, Tc]);
            obj.solution{2} = chebfun(my(tList), [0, Tc]);
            obj.solution{3} = chebfun(mz(tList), [0, Tc]);
        end
        

    end
    
    methods (Access = private)
        [mx, my, mx_list, my_list] = mxmy(obj, bx, by, bz, t, x, zc, varargin)
        [ac, as, fc, fs] = amp_phase(obj, x, z0, zc, p, nmax)
        [mx, my] = m_p(obj, bx, by, t, x, z0, zc, p, nmax)        
        [Wp, Wm, Qp, Qm] = WQ(obj, x, z0, zc, p, nmax)
        [Wp, Wm, Qp, Qm] = WQ_n(obj, x, z0, zc, n, p)        
    end
end

